5 Diagnósticos de Residuos

5.1 Pretratamientos de la Información

Cargar Recursos Externos

source("scripts/recursos.R")
source("scripts/funciones.R")

Definición de variables generales

## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <- 
  "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Definición de variables específicas

# List of crimes/delitos and cities/ciudades
crimes = c(
  "Minor property offenses",
  "Robberies",
  "Burglaries",
  "Drunkennes, damages and disorders",
  "Family violence and aggression",
  "Injuries, drugs and weapons",
  "High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")

#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]

Lectura de Datos Espaciales

Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.

city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))

El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame.

Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)

# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)

Transformaciones a los datos.

Se utilizó funciónes de transformacion general mencionadas en Section 6.2 que corresponden a las siguientes:

Transformación Cúbica
Transformación Logarítmica
Cáculo del Lag Espacial

cityp <- cityp %>%
    st_as_sf(crs = 4326) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_cub, .names = "cub{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_log1p, .names = "log{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = lag_spatial, nb, .names = "lag{.col}")) %>%
    as("Spatial")

zona	id	x	y	pflot	abus	amen	asalto	denspob	educ	tamhog	constr	ofcom	vulzon	mobil	lpflot	neduc	ldenspob	lofcom	lconstr	cubasalto	cubamen	cubabus	logasalto	logamen	logabus	laghurto	lagrobo	lagasalto	lagamen	lagabus	lagcrim
4101161005	22	-71.239	-29.862	1.000	0	0	0	144.756	10.867	3.170	4328.618	0.000	0.367	0.079	0.178	0.477	0.767	0.000	0.840	0	0	0	0.000	0.000	0.000	0.083	0.083	0.083	0.417	0.250	0.083
4101161005	32	-71.243	-29.863	1.000	0	0	0	227.772	10.716	3.393	3830.009	0.000	0.367	0.079	0.178	0.466	0.846	0.000	0.829	0	0	0	0.000	0.000	0.000	0.000	0.083	0.083	0.917	0.583	0.083
4101161005	35	-71.240	-29.863	1.034	0	0	1	171.133	11.801	3.172	4328.618	0.000	0.367	0.079	0.182	0.550	0.796	0.000	0.840	1	0	0	0.693	0.000	0.000	0.000	0.083	0.000	0.500	0.250	0.083
4101161005	36	-71.239	-29.863	1.000	0	0	0	154.603	10.870	3.135	4328.618	0.000	0.367	0.079	0.178	0.478	0.779	0.000	0.840	0	0	0	0.000	0.000	0.000	0.083	0.083	0.083	0.417	0.250	0.083
4101161005	37	-71.238	-29.863	1.025	0	0	0	155.145	10.858	3.148	3230.309	0.000	0.367	0.079	0.181	0.477	0.779	0.000	0.815	0	0	0	0.000	0.000	0.000	0.083	0.083	0.083	0.000	0.250	0.083
4101161001	48	-71.244	-29.864	1.125	0	0	0	168.179	11.958	3.336	2384.222	14.682	0.292	0.122	0.193	0.562	0.793	0.574	0.788	0	0	0	0.000	0.000	0.000	0.000	0.000	0.250	0.500	0.417	0.083
4101161005	49	-71.243	-29.864	1.071	0	0	0	212.743	10.978	3.415	3795.759	0.000	0.367	0.079	0.187	0.486	0.834	0.000	0.828	0	0	0	0.000	0.000	0.000	0.000	0.083	0.250	0.917	0.583	0.083
4101161005	50	-71.242	-29.864	1.000	1	1	0	220.188	10.705	3.409	4026.613	0.000	0.367	0.079	0.178	0.465	0.840	0.000	0.834	0	1	1	0.000	0.693	0.693	0.000	0.083	0.000	1.000	0.417	0.083
4101161005	51	-71.241	-29.864	1.021	1	1	0	243.501	10.120	3.337	4328.618	0.000	0.367	0.079	0.180	0.419	0.857	0.000	0.840	0	1	1	0.000	0.693	0.693	0.000	0.083	0.083	0.667	0.250	0.000
4101161005	52	-71.240	-29.864	1.022	0	0	0	320.023	10.720	3.324	4328.618	0.000	0.367	0.079	0.180	0.466	0.905	0.000	0.840	0	0	0	0.000	0.000	0.000	0.000	0.083	0.167	0.583	0.250	0.000

5.2 Evaluación De Residuos

5.2.1 Modelos MLB

Lectura del Modelo.

Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string

Warning in wkt(obj): CRS object has no comment

model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))

5.2.1.1 Extracción de Residuos

cityp$res_fit_mlb <- residuals(model_mlb)
cityp$fitted_fit_mlb <- fitted(model_mlb)
cityp$sd_breaks_mlb <- scale(cityp$res_fit_mlb)[,1]

breaks_qt <- classInt::classIntervals(cityp$sd_breaks_mlb, 
                                      n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)

city_sf <- cityp %>% st_as_sf(crs = 32719) 
residuos_mv <- mapview(city_sf,zcol="res_fit_mlb", at =my_breaks, cex=2) 
residuos_mv

5.2.1.2 Global Moran

Global Moran a Variable Depediente

spdep::moran.test(x = cityp$laghurto,  listw = nb)


    Moran I test under randomisation

data:  cityp$laghurto  
weights: nb    

Moran I statistic standard deviate = 174.4, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.7640052224     -0.0001215805      0.0000191963

# moran.plot(x =cityp$laghurto,  listw = nb, labels=as.character(cityp$id))

Global Moran de Residuos

spdep::moran.test(x = cityp$res_fit_mlb,  listw = nb)


    Moran I test under randomisation

data:  cityp$res_fit_mlb  
weights: nb    

Moran I statistic standard deviate = 19.256, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
    0.08478830236    -0.00012158055     0.00001944476

# moran.plot(x =cityp$res_fit_mlb,  listw = nb, labels=as.character(cityp$id))

5.2.1.3 Local Moran de Residuos

# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_mlb, nb)
cityp$res_scaled = as.numeric(scale(cityp$res_fit_mlb))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])

# Umbral de significancia estadistica
pval=0.05

# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)


  HH   HL   LH   LL   NS 
 284   33  295  331 7283

city_sf <- cityp %>% st_as_sf(crs = 32719) %>% 
  filter(clusterM != "NS")

residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran

5.2.2 Modelos SAR

Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string

Warning in wkt(obj): CRS object has no comment

model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))

5.2.2.1 Extracción de Residuos SAR

cityp$res_fit_sar <- residuals(model_sar)
cityp$fitted_fit_sar <- fitted(model_sar)

This method assumes the response is known - see manual page

cityp$sd_breaks_sar <- scale(cityp$res_fit_sar)[,1]

breaks_qt <- classInt::classIntervals(cityp$sd_breaks_sar, 
                                      n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)

city_sf <- cityp %>% st_as_sf(crs = 32719) 

residuos_mv <- mapview(city_sf,zcol="res_fit_sar", at =my_breaks, cex=2) 
  
residuos_mv

5.2.2.2 Global Moran

Global Moran a Variable Depediente

spdep::moran.test(x = cityp$cubhurto,  listw = nb)


    Moran I test under randomisation

data:  cityp$cubhurto  
weights: nb    

Moran I statistic standard deviate = 79.737, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
    0.35153153375    -0.00012158055     0.00001944977

# moran.plot(x =cityp$laghurto,  listw = nb, labels=as.character(cityp$id))

Global Moran de Residuos

spdep::moran.test(x = cityp$res_fit_sar,  listw = nb)


    Moran I test under randomisation

data:  cityp$res_fit_sar  
weights: nb    

Moran I statistic standard deviate = -1.2062, p-value = 0.8861
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
   -0.00544188551    -0.00012158055     0.00001945508

# moran.plot(x =cityp$res_fit_mlb,  listw = nb, labels=as.character(cityp$id))

5.2.2.3 Local Moran de Residuos

# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_sar, nb)

cityp$res_scaled = as.numeric(scale(cityp$res_fit_sar))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])

# Umbral de significancia estadistica
pval=0.05

# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)


  HH   HL   LH   LL   NS 
 214   19  156   32 7805

city_sf <- cityp %>% st_as_sf(crs = 32719) %>% 
  filter(clusterM != "NS")

residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran